Problem Statement¶
WHO wants a data driven approach which could help in suggesting a country which area should be given importance in order to efficiently improve the life expectancy of its population. Life_Expectancy_Data.csv file contains the data
Solution¶
Building a Prediction engine which predicts the life expectancy based on various features like status of the country, GDP, Alcohol consumption, adult mortality rate etc. To build a prediction engine we will use linear regression model.
Purpose: To perform regression analysis to study life expectancy in different countries
Linear Regression for the above business problem.¶
- Loading the Data.
- Understanding the Data
- Data Preprocessing
- Exploratory Data Analysis
- Model Building
- Model Diagnostics
- Predictions and Evaluations
Importing Libraries- Setting up the environment¶
- Libraries Numpy and Panda are used for data manipulation and data loading.
- Libraries Seaborn and Matplotlib are used for data visualization.
- Libraries Sklearn and Statsmodel for data preprocessing and model building.
In [1]:
import numpy as np
import pandas as pd
pd.set_option('display.max_rows',800)
pd.set_option('display.max_columns',500)
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline
# importing all libraries and dependencies for machine learning
from sklearn import preprocessing
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from sklearn.feature_selection import RFE
from sklearn.linear_model import LinearRegression
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.metrics import mean_absolute_error, mean_squared_error,r2_score
import random
1. Loading the data¶
In [2]:
df = pd.read_csv("~/data/Life_Expectancy_Data.csv")
2.Understanding the data¶
In [3]:
df.info()
In [4]:
df.describe()
Out[4]:
In [5]:
df.head(5)
Out[5]:
In [6]:
num_col = df.select_dtypes(include=np.number).columns
print("Numerica columns are as follows.\n",num_col)
cat_col = df.select_dtypes(exclude=np.number).columns
print("Categorical columns are as follows.\n",cat_col)
3. Data Pre-Processing¶
In [7]:
#Removing the extra space from column names.
#Strip will remove spaces before and after the column names,if any.
#Using Lambda helps scaling up to all the columns.
df = df.rename(columns = lambda x : x.strip())
In [8]:
#converting categorical features to numerical features for the machine to understand.
from sklearn import preprocessing
#Label encoder understands labels
label_encoder = preprocessing.LabelEncoder()
#Status "Developing" will be 1 and "Developed" will be 0
df['Status'] = label_encoder.fit_transform(df['Status'])
df.head()
Out[8]:
In [9]:
#Knowing Missing values
print(df.isna().sum())
print(df.shape)
In [10]:
# Treating the NA values
# Replacing the NA value with Mean of the column
for i in df.columns.drop('Country'):
df[i].fillna(df[i].mean(),inplace = True)
4. Exploratory Data Analysis¶
In [11]:
#Checking the distribution of Y variable.
plt.figure(figsize=(8,8),dpi=80)
sns.boxplot(df['Life expectancy'])
plt.title('Life Exp Box plot')
plt.show()
In [12]:
plt.figure(figsize=(8,8))
plt.title('Life Exp Distribution plot')
sns.distplot(df['Life expectancy'])
Out[12]:
Summary : The Y variable has few outliers and is almost linearlly distributed. Hence the Assumption of the linear regression holds true.
In [13]:
num_col = df.select_dtypes(include = np.number).columns
cat_col = df.select_dtypes(exclude = np.number).columns
print('Numerical columns are ',num_col)
print('categorical columns are ',cat_col)
In [14]:
#checking the multicolleaniarity of features by checking the correlation matrix
plt.figure(figsize=(16,16))
p = sns.heatmap(df[num_col].corr(), annot = True, cmap = 'RdYlGn', center = 0)
In [15]:
# To know the relation between different features
ax = sns.pairplot(df[num_col])
5. Model Building¶
In [16]:
# Train test Split (70% data - Train, 30% data - Test)
X= df.drop(columns=['Life expectancy','Country'])
y=df[['Life expectancy']]
X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.3, random_state=1234)
Approach 1: Adding 1 variable at a time.¶
Building model with 1 variable¶
In [17]:
X_train1 = X_train['Income composition of resources']
In [18]:
# Adding a constant
X_train1 = sm.add_constant(X_train1)
# Creating a OLS model
model_1 = sm.OLS(y_train, X_train1).fit()
In [19]:
model_1.params
Out[19]:
In [20]:
print(model_1.summary())
Building model with 2 variables¶
In [21]:
#Adding Schooling feature to the regression model
X_train2 = X_train[['Income composition of resources','Schooling']]
In [22]:
#Adding constant
X_train2 = sm.add_constant(X_train2)
In [23]:
model_2 = sm.OLS(y_train, X_train2).fit()
In [24]:
model_2.params
Out[24]:
In [25]:
print(model_2.summary())
Building with 3 variables¶
In [26]:
# adding 3rd feature in regression model
X_train3 = X_train[['Income composition of resources', 'Schooling', 'Adult Mortality']]
In [27]:
#adding a constant
X_train3 = sm.add_constant(X_train3)
In [28]:
#creating fitted model
model_3 = sm.OLS(y_train, X_train3).fit()
In [29]:
#checking the parameters
model_3.params
Out[29]:
In [30]:
#Printing the model summary
print(model_3.summary())
Approach 2: Solving by Recursive Feature elimination(RFE) and eliminating by using P-Value and VIF¶
In [31]:
#Running RFE with important columns(count 15 imp. features)
lm = LinearRegression()
lm.fit(X_train,y_train)
rfe = RFE(lm, 15)
rfe = rfe .fit(X_train, y_train)
In [32]:
list(zip(X_train.columns, rfe.support_,rfe.ranking_))
Out[32]:
In [33]:
#Selecting important features basis support
imp_columns = X_train.columns[rfe.support_]
In [34]:
print(imp_columns)
In [35]:
#creating X_train dataframe with RFE selected variables
X_train_rfe = X_train[imp_columns]
After passing selected columns by RFE, manually evaluate each models p-value and VIF value. Until the acceptable range for p-values and VIF is found, drop the variables one at a time basis below criteria.
- High VIF High p-value : Drop the variable
- Low VIF High p-value : Drop the variable with high p-value first
- Low VIF Low p-value : Accept the variable
In [36]:
random.seed(0)
#Add constant
X_train_rfec = sm.add_constant(X_train_rfe)
#model with RFE features
lm_rfe = sm.OLS(y_train, X_train_rfec).fit()
print(lm_rfe.summary())
In [37]:
vif = pd.DataFrame()
vif['Features'] = X_train_rfe.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe.values, i) for i in range(X_train_rfe.shape[1])]
vif['VIF'] = round(vif['VIF'], 2)
vif = vif.sort_values(by = 'VIF', ascending = False)
vif
Out[37]:
Dropping 'thinness 5 - 9 yrs' feature from training dataset as it has high p value.
In [38]:
X_train_rfe1 = X_train_rfe.drop(['thinness 5-9 years'],1,)
#add constant
X_train_rfe1b = sm.add_constant(X_train_rfe1)
#building model
lm_rfe1 = sm.OLS(y_train, X_train_rfe1b).fit()
#summary
print(lm_rfe1.summary())
In [39]:
#creating dataframe that contains feature names and its VIF
vif = pd.DataFrame()
vif['Features'] = X_train_rfe1.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe1.values,i) for i in range(X_train_rfe1.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif= vif.sort_values(by ='VIF', ascending = False)
vif
Out[39]:
Under five deaths has high VIF hence dropping that feature
In [40]:
X_train_rfe2 = X_train_rfe1.drop('under-five deaths',1,)
#adding constant
X_train_rfe2c = sm.add_constant(X_train_rfe2)
#building model
lm_rfe2 = sm.OLS(y_train, X_train_rfe2c).fit()
#summary
print(lm_rfe2.summary())
In [41]:
vif = pd.DataFrame()
vif['Features'] = X_train_rfe2.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe2.values,i) for i in range(X_train_rfe2.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif = vif.sort_values(by = 'VIF', ascending = False)
vif
Out[41]:
Since variable Alcohol has a very high P value eliminating that feature from the training dataset
In [42]:
#Dropping insignificant variable Alcohol
X_train_rfe3 = X_train_rfe2.drop('Alcohol',1,)
#adding constant
X_train_rfe3c = sm.add_constant(X_train_rfe3)
#model building
lm_rfe3 = sm.OLS(y_train, X_train_rfe3c).fit()
#summary
print(lm_rfe3.summary())
In [43]:
#VIF for all features
vif = pd.DataFrame()
vif['Features'] = X_train_rfe3.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe3.values,i) for i in range(X_train_rfe3.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif = vif.sort_values(by = 'VIF', ascending = False)
vif
Out[43]:
Dropping 'Schooling' Feature from training data set since it has high VIF
In [44]:
X_train_rfe4 = X_train_rfe3.drop('Schooling',1,)
#add a constant
X_train_rfe4c = sm.add_constant(X_train_rfe4)
#Model building
lm_rfe4 = sm.OLS(y_train, X_train_rfe4c).fit()
#Summary
print(lm_rfe4.summary())
In [45]:
vif = pd.DataFrame()
vif['Features'] = X_train_rfe4.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe4.values,i) for i in range(X_train_rfe4.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif = vif.sort_values(by = 'VIF', ascending = False)
vif
Out[45]:
'Diphtheria' has a high VIF value hence dropping the feature from the training set.
In [46]:
X_train_rfe5 = X_train_rfe4.drop('Diphtheria',1,)
#adding constant
X_train_rfe5c = sm.add_constant(X_train_rfe5)
#building model
lm_rfe5 = sm.OLS(y_train, X_train_rfe5c).fit()
#summary
print(lm_rfe5.summary())
In [47]:
vif = pd.DataFrame()
vif['Features'] = X_train_rfe5.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe5.values,i) for i in range(X_train_rfe5.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif = vif.sort_values(by ='VIF', ascending = False)
vif
Out[47]:
'Polio' has a high VIF value hence dropping the feature from the training set.
In [48]:
X_train_rfe6 = X_train_rfe5.drop('Polio',1,)
#adding constant
X_train_rfe6c = sm.add_constant(X_train_rfe6)
#building model
lm_rfe6 = sm.OLS(y_train, X_train_rfe6c).fit()
#summary
print(lm_rfe6.summary())
In [49]:
vif = pd.DataFrame()
vif['Features'] = X_train_rfe6.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe6.values,i) for i in range(X_train_rfe6.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif = vif.sort_values(by ='VIF', ascending = False)
vif
Out[49]:
'Hepatitis B' has a high VIF value hence dropping the feature from the training set.
In [50]:
X_train_rfe7 = X_train_rfe6.drop('Hepatitis B',1,)
#adding constant
X_train_rfe7c = sm.add_constant(X_train_rfe7)
#building model
lm_rfe7 = sm.OLS(y_train, X_train_rfe7c).fit()
#summary
print(lm_rfe7.summary())
In [51]:
vif = pd.DataFrame()
vif['Features'] = X_train_rfe7.columns
vif['VIF'] = [variance_inflation_factor(X_train_rfe7.values,i) for i in range(X_train_rfe7.shape[1])]
vif['VIF'] = round(vif['VIF'],2)
vif = vif.sort_values(by ='VIF', ascending = False)
vif
Out[51]:
Approach 3: Stepwise Regression¶
In [52]:
## By David Dale https://datascience.stackexchange.com/users/24162/david-dale
def stepwise_selection(X, y,
initial_list=[],
threshold_in=0.05,
threshold_out = 0.10,
verbose=True):
""" Perform a forward-backward feature selection
based on p-value from statsmodels.api.OLS
Arguments:
X - pandas.DataFrame with candidate features
y - list-like with the target
initial_list - list of features to start with (column names of X)
threshold_in - include a feature if its p-value < threshold_in
threshold_out - exclude a feature if its p-value > threshold_out
verbose - whether to print the sequence of inclusions and exclusions
Returns: list of selected features
Always set threshold_in < threshold_out to avoid infinite looping.
See https://en.wikipedia.org/wiki/Stepwise_regression for the details
"""
included = list(initial_list)
while True:
changed=False
# forward step
excluded = list(set(X.columns)-set(included))
new_pval = pd.Series(index=excluded)
for new_column in excluded:
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included+[new_column]]))).fit()
new_pval[new_column] = model.pvalues[new_column]
best_pval = new_pval.min()
if best_pval < threshold_in:
best_feature = new_pval.idxmin()
included.append(best_feature)
changed=True
if verbose:
print('Add {:30} with p-value {:.6}'.format(best_feature, best_pval))
# backward step
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
# use all coefs except intercept
pvalues = model.pvalues.iloc[1:]
worst_pval = pvalues.max() # null if pvalues is empty
if worst_pval > threshold_out:
changed=True
worst_feature = pvalues.argmax()
included.remove(worst_feature)
if verbose:
print('Drop {:30} with p-value {:.6}'.format(worst_feature, worst_pval))
if not changed:
break
return included
result = stepwise_selection(X_train, y_train)
print('resulting features:')
print(result)
In [53]:
X_train_stepwise = X_train[['Schooling', 'Adult Mortality', 'HIV/AIDS', 'Diphtheria', 'BMI', 'Income composition of resources', 'Status', 'percentage expenditure', 'Polio', 'Measles', 'Hepatitis B', 'under-five deaths', 'infant deaths', 'thinness 1-19 years']]
# Adding constant
X_train_stepwise = sm.add_constant(X_train_stepwise)
# Building a model
lm_stepwise = sm.OLS(y_train, X_train_stepwise).fit()
# Summary
print(lm_stepwise.summary())
Model Prediction and Evaluation¶
In [56]:
# Prediction with Training set
X_test_stepwise = X_test[['Schooling', 'Adult Mortality', 'HIV/AIDS', 'Diphtheria', 'BMI', 'Income composition of resources', 'Status', 'percentage expenditure', 'Polio', 'Measles', 'Hepatitis B', 'under-five deaths', 'infant deaths', 'thinness 1-19 years']]
X_test_stepwise = sm.add_constant(X_test_stepwise)
actual = y_test['Life expectancy']
prediction = lm_stepwise.predict(X_test_stepwise)
In [57]:
# Model Evaluation: Calculating Mean Squared Errors
model_mse = mean_squared_error(prediction, actual)
print(model_mse)
In [58]:
# Model Evaluation: Calculating Mean Absolute percentage Errors
def mape(y_true, y_pred):
y_true, y_pred = np.array(y_true), np.array(y_pred)
return np.mean(np.abs((y_true - y_pred) / y_true)) *100
In [59]:
mape(actual,prediction)
Out[59]:
In [60]:
# Linearity check
sns.scatterplot(y_test['Life expectancy'],prediction)
plt.title('Linearity Check')
plt.xlabel('Actual value')
plt.ylabel('Predicted value')
Out[60]:
In [61]:
# Histogram to analyse error terms
fig = plt.figure()
sns.distplot((y_test['Life expectancy'] - prediction), bins = 20)
fig.suptitle('Error Term Analysis', fontsize = 20)
plt.xlabel('Errors', fontsize = 18)
Out[61]:
In [ ]: